Anomalous low doping phase of the Hubbard model 



C. Grober, R. Eder and W. Hanke 
Institut fiir Theoretische Physik, Universitdt Wiirzburg, Am Hubland, 97074 Wiirzburg, Germany 

(February 1, 2008) 

We present results of a systematic Quantum-Monte-Carlo study for the single-band Hubbard 
model. Thereby we evaluated single-particle spectra (PES & IPES), two-particle spectra (spin & 
density correlation functions) , and the dynamical correlation function of suitably defined diagnostic 
operators, all as a function of temperature and hole doping. The results allow to identify differ- 
ent physical regimes. Near half-filling we find an anomalous 'Hubbard-I phase', where the band 
structure is, up to some minor modifications, consistent with the Hubbard-I predictions. At lower 
temperatures, where the spin response becomes sharp, additional dispersionless 'bands' emerge due 
to the dressing of electrons/holes with spin excitatons. We present a simple phenomenological fit 
which reproduces the band structure of the insulator quantitatively. The Fermi surface volume in 
the low doping phase, as derived from the single-particle spectral function, is not consistent with the 
Luttinger theorem, but qualitatively in agreement with the predictions of the Hubbard-I approxi- 
mation. The anomalous phase extends up to a hole concentration of ~ 15%, i.e. the underdoped 
region in the phase diagram of high-Tc superconductors. We also investigate the nature of the 
magnetic ordering transition in the single particle spectra. We show that the transition to an SDW- 
like band structure is not accomplished by the formation of any resolvable 'precursor bands', but 
rather by a (spectroscopically invisible) band of spin 3/2 quasiparticles. We discuss implications for 
the 'remnant Fermi surface' in insulating cuprate compounds and the shadow bands in the doped 
materials. 
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I. INTRODUCTION 

The 1-band Hubbard model on a two-dimensional 
square lattice has the Hamiltonian 

i 

(1) 



Here, c] ^ i'^i a) creates (annihilates) an electron with 
spin cr in a Wannier orbital centered at lattice site i. The 
particle density at each site is given by rii^^ — c\^c^^. 
The first sum for the kinetic energy is restricted to in- 
clude only the hopping matrix element t between next- 
nearest neighbor sites (i, j). Periodic boundary condi- 
tions are used throughout the following work. The sec- 
ond sum describes for [/ > an on-site Coulomb re- 
pulsion between particles of opposite spin that share the 
same lattice site. In the present work we restrict our- 
selves to U = 8.0 1. The chemical potential /i in the 
third sum controls the occupation of the finite lattice in 
the finite-temperature grand canonical Quantum-Monte- 
Carlo (QMC) simulation we preformed. At half-filling, 
particle-hole symmetry of the kinetic and [/-term im- 
plies /i = 0. The analytic continuation of the dynamic 
imaginary-times QMC data to the real frequency axis is 
performed with state-of-the-art Maximum-Entropy (ME) 
techniques. For exhausting discussions concerning the 
QMC and ME methods we refer the reader to pl-ll. 



The 1-band Hubbard model exhibits several energy 
scales: in the repulsive case the high-energy scale U 
is important in determining the insulating gap at half- 
filling, (n) = 1.0. An important low-energy scale is set 
by the exchange interaction J = At^ /U: in second order 
perturbation theory two particles with different spins on 
neighboring lattice sites can exchange via a virtual dou- 
ble occupation. This process is the source for the strong 
antiferromagnetic (AF) correlations found near and at 
half-filling. 

With the exception of some known symmetry properties 
like invariance under global spin-rotation, whose gen- 
erators form the S'C/(2)-algebra, and invariance under 
[/(l)-transformation, i.e., charge-conservation, as well as 
the particle-hole transformation, no rigorous results are 
known for the 1-band Hubbard model in two-dimensions 
Q. The Mermin- Wagner theorem prevents a long- 
range ordered state in a two-dimensional system for fi- 
nite temperatures, but it is commonly believed that the 
ground state of the spin-1/2 Heisenberg antiferromagnet, 
i.e. the large-C/ limit of the repulsive half-filled Hubbard 
model 1^ , shows long-range Neel order in two dimensions. 
Neel order results also in the weak-coupling limit [0,D, 
where the gap A is due to a spin-density-wave (SDW) 
instability related to perfect nesting. 
It might seem that due to the Mermin- Wagner theorem 
the physics of the ordered phase is out of reach for our 
numerical technique, which is limited to finite temper- 
atures. However, one may assume that the system 'ef- 
fectively orders' as soon as the spin-correlation length 
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becomes comparable to the system size. Due to the pe- 
riodic boundary conditions the spin-correlation function 
x{r) — ■ ^i+r) is a periodic function of r 

and if the value of this function at the maximum value 
|r| = ^/2L (with L = cluster size) is still appreciable, 
we may expect that the system is 'effectively ordered'. 
A rough measure would be the spin-correlation length ^, 
obtained by fitting xi'"') to the form a ■ \r\'' ■ e"''''/^. Since 
the infinite system has the Neel temperature T ~ 0, the 
spin-correlation length diverges as T — > and we may ex- 
pect that in a finite system C becomes comparable to the 
cluster size L at a finite temperature which depends on 
the lattice size. Below this temperature we then expect 
that the system resembles the ordered phase, although 
the spin-rotation symmetry persists even in this case. In 
that sense, the finite size of the system creates an arti- 
ficial 'Neel temperature' which, however, depends on L 
and U/t etc and has no real counterpart in the infinite 
system. This has to be kept in mind when discussing the 
results. 

We proceed by discussing some known approximations 
to the Hubbard model. The 'classical' approximation to 
the Hubbard model is the so-called Hubbard-I approxi- 
mation |10|. Its essence is the splitting of the electron 
annihilation operator into the two 'eigenoperators' of the 
interaction part: 



of these effective particles, and an additional energy of 
formation of U for the double occupancy. The matrix 
elements for the propagation are reduced by a factor of 
1/2, because in an uncorrelated spin-background there 
is a probability of 1/2 for the spin on a nearest neigh- 
bor to have the proper direction to allow for the hopping 
of an electron. Solving (^) by Fourier and Bogoliubov 
transform: 



l2M,a = -Vkdk.a + ^khl^g, 



(5) 



yields the standard dispersion relation, which consists of 
the upper and lower Hubbard bands: 



E. 



Hub-I 



(k) 



efc ± 



(6) 



Here efc denotes the free tight-binding dispersion efc = 
—2t{cos{kx) + cos{ky)). Using (||) we also obtain the cor- 
rect spectral weights of the two Hubbard bands: 



(k) = -(ufc ± VkY 



1 / i , <^fe 



^ 1± 



(7) 
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whence 
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The physical content of the Hubbard-I approximation, 
which neglects any spin-correlations, becomes clear by 
realizing pT] that the equations of motion in this 
approximation are completely equivalent to an 'effec- 
tive Hamiltonian' for double occupancy-like particles 
d] ^ = (l/\/2) cj ^Ui^g and hole-hke particles hi ^ = 
(1/V2) c,,,(l-n-s): 



E idlhl + h.c.) 



U 



(4) 



This Hamiltonian contains terms which describe the pair 
creation of a hole and a double occupancy on nearest 
neighbors terms which describe the propagation 



As noted above, the key assumption in the Hubbard I 
approximation is the neglect of spin-correlations. There- 
fore, we expect that this approximation will become in- 
accurate as soon as spin-correlations become sufficiently 
strong so as to appreciably influence the propagation 
of holes and double occupancies. This effect will be 
strongest in the Neel ordered phase believed to be re- 
alized in the ground state. If we choose the 'spin- 
background', in which the holes and double occupancies 
propagate, to be the Neel state, the double occupancy 
c?|-|. with spin | can exist only on the J, sublattice and 

vice versa. Similarly, a hole can be created only on 
the I sublattice and vice versa. We, thus, expect that we 
have to modifiy the Hubbard-I Hamiltonian (Q) into: 



H 



E {dl-^ht.^ + h.c.) 

ieA,jeN{i) 



U 



(8) 



where A [B] denotes the | (t) sublattice and N{i) de- 
notes the set of nearest neighbors of site i. Note that 
the d\^dj,a and h\ ^hj^^ propagation terms drop out (be- 
cause of the A and B sublattices), and the matrix element 
for pair creation are t rather than t/2 - this takes into 
account the fact that the spins on neighboring sites are 
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antiparallel with probability 1 rather than 1/2 as was the 
case in the paramagnetic state. Fourier and Bogohubov 
transformation of (g) yields the dispersion 



and using 



(9) 



(10) 



(where k is within the AF Brillouin zone), we find the 
spectral weight of these bands: 



Z± = ^{uk ± Vk f 




(11) 



This is precisely what is obtained from the SDW mean- 
field treatment of the Hubbard model by setting the 
staggered magnetization m to a value of 1 (which is a 
good approximation in the limit of large U). In general, 
the SDW mean-field approximation models the AF Neel 
state by assuming {ni,a) — 5(1 + crme**^^') with AF 
nesting vector Q = (tt, tt) and staggered magnetization 



dispersion 



)|| 



This results in the two-band 



i?|^^(k) = ±^4-|-A2, 
and the spectral weight 

^r^(k) = i(l±Wi?r^(k)) 



(12) 



(13) 



The gap parameter is A — Um/2 and the staggered mag- 
netization m is determined self-consistently from the fol- 
lowing equation: 



m 



occupied 



Um 



N ^ V(efe-e(k + Q))2+m2[/2 



(14) 



The solution of this self-consistency equation yields the 
value A = 3.56 i at C/ = 8.0 1, the value used in the 
present work. If we set m = 1 on the other hand, as 
would be apropriate for V jt 00, we obviously recover 
the results from our Hubbard- I-like Hamiltonian (H). 
In the two limiting cases of no spin-correlations and of 
perfect Neel order we can thus treat the Hubbard model 
in a quite analogous fashion and, as will be seen be- 
low, the Hubbard-I results are indeed a good approx- 
imation to the actual spectral function in the limit of 
high temperature. The main problem then is how to de- 
scribe the effect of spin-fluctuations and how to manage 



the crossover from the completely disordered to the Neel 
ordered phase. Below, we will adress this crossover by 
QMC simulations. 

Thereby we want to take advantage of the possibility to 
calculate the spectra of specifically designed 'diagnostic 
operators'. The first one of these is the 'shadow' operator. 



(15) 



We note first of all, that the Fourier transform of this 
operator, Ck^a- has precisely the same quantum numbers 
as the ordinary electron operator Ck.a'- momentum — fc, 
total spin 1/2, z-spin cr, and identical point group sym- 
metry at high symmetry momenta. It follows that the 
poles in its dynamical correlation functions 

Aa(fc,o;)=5]|(vI',|5fc^,|vI/^,)|2 

originate from exactly the same final states j^'^) as those 
of the photoemission spectrum. It can happen only ac- 
cidentally that a given state v has an exactly vanish- 
ing weight in one of the correlation functions, but not 
the others. In this case, however, any arbitrarily small 
perturbation will remove the accidental vanishing of the 
peak. The only thing that can and will be different in the 
spectra of the diagnostic operator, are the weights of the 
peaks, |(^'i/|cfc,o-|^o)P- In fact, comparison of equations 
( |TT| ) and (|^) shows immediately that both for the Hub- 
bard I approximation and for the SDW-approximation 
use of the operator Ck.a instead of Ck.a exchanges the 
weights of the two Hubbard/SDW bands. It follows that 
band portions which have a small spectral weight in the 
spectra of the ordinary electron operator will aquire a 
large spectral weight in the spectrum of Ck.a and vice 
versa. Therefore, this diagnostic operator is a useful tool 
to map out the 'shadowy' parts of the spectra in the outer 
parts of the Brillouin zone. An additional benefit is that 
since the ME technique resolves peaks with large spec- 
tral weight more reliably than those with small weight, 
we can get more precise information about the dispersion 
of these bands with weak intensity. The second diagnos- 
tic operator that we will be using is the 'spin-1/2 string' 
operator 



( + CiX'S'j ), 

3eN(i) 



(17) 



where the sum on the r.h.s. runs over the nearest neigh- 
bors j of site i. This operator is the Clebsch-Gordan con- 
traction of the adjoint rank-1 spinor Ci^^ and the spin-1 
vector operator Sj into yet another adjoint rank-1 spinor. 
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It describes the annihilation of a 'dressed' electron, i.e., 
an electron with a spin-excitation on a nearest neighbor. 
Again, since the Fourier transforms of the diagnostic op- 
erator, bk^a-i agrees with the electron annihilation opera- 
tor Ck,a in all conceivable quantum numbers (momentum, 
total spin, z-component of the spin, point group opera- 
tions in the case of high symmetry momenta), it obeys 
precisely the same selection rules, whence it is just as 
good as the c-operator itself to map out the band struc- 
ture. 

One further point, which is important from the techni- 
cal point of view, is the following: under the particle- 
hole transformation Q o- — > e^^'-^'c] we have c,- 



and b. 



giQHi^t^^ This implies that at 



half- filling the spectra of the shadow and spin- 1/2 string 
operator obey the same particle-hole-symmetry as those 
of the ordinary electron operator, i.e. 



A{k,uj) = A{k + Q, -to) 



(18) 



In our MaxEnt program particlc-holc-symmctry is not 
implemented as an additional constraint, in other words: 
the MaxEnt procedure does not 'know' about this addi- 



tional symmetry. The degree to which (18) is obeyed in 
the final spectra thus gives a good check for the accuracy 
of the reconstructed spectra. This is of particular impor- 
tance in the case of the spin-1/2 string operator because 
the Wick-contraction of this operator on any given time 
slice produces a total of « 80 products of noninteracting 
Green's functions. The computation is, therefore, much 
more prone to inaccuracies so that an additional check is 
desirable. 

The present work is organized as follows: first, we com- 
pare the temperature dependent dynamic single-particle 
properties of the Hubbard model with the predicitions of 
the mean-field SDW and Hubbard-I approximations. In 
addition, we consider the temperature dependent two- 
particle excitations. Then, we will use our first diag- 
nostic operator, the shadow operator Ck^a, to show the 
existence of a total of four bands in the photoemission 
spectrum and to shed light onto the temperature depen- 
dent crossover from the SDW to the Hubbard-I regime. 
Then we will investigate the 4-band structure in more 
detail: after a phenomenological fit, which is able to pro- 
duce a total of four bands, we will consider the string 
picture which naturally leads to our second diagnostic 
operator, the spin-1/2 string operator bk,a- This op- 
erator will be used to ultimately reveal the underlying 
mechanisms behind the generation of the 4-band struc- 
ture, namely the dressing of the photoholes by clouds of 
AF spin-excitations. To resolve the 'AF mirror image' of 
the narrow quasiparticle spectral weight features between 
k = (0, 0) and k = (tt, 0) around momentum k = (tt, tt) 
we introduce also a spin-3/2 string operator. Finally, we 
will concentrate on the doped regime, thereby showing 
the violation of the Luttinger theorem near half-filling, 



and discuss the hole concentration range in which these 
dressing effects dominante the low-energy physics. 



II. TEMPERATURE DEPENDENT DYNAMICS 
AT HALF-FILLING 

We start in figure ^ with the discussion of the angle- 
resolved single-particle spectral function A(k, oj) for U = 
8.0 1 and various temperatures in the range from T = 
4.00 1 to T = 0.10 1. In the left column of the figure 
the spectral functions are shown as 'grey-scale' plots ver- 
sus momentum k and energy uj/t with dark (light) ar- 
eas corresponding to large (small) spectral weight. The 
same spectra are shown also in the right column, but 
now as line-plots at each momentum k. The QMC data 
in the figure are compared to the renormalized results 
of the Hubbard-I approximation at high and medium 
temperatures and with the renormalized results of the 
mean-field SDW approximation at the lowest tempera- 
ture, T = 0.10 1. 'Renormalized' means here that we have 
readjusted the parameters U and t in ^ and (p^so as 
to obtain an optimal fit to the 'bands' of high spectral 
weight in the spectra. These approximate dispersions are 
plotted as solid lines in the left column, while their spec- 
tra are shown in the right column as line-plots at each 
momentum k. Thereby we have assumed a Lorentzian 
lineshape with a suitably chosen temperature dependent 
width (the Hubbard-I approximation does not provide 
any information about linewidths). 

Starting at the highest calculated temperature, T — 
4.00 1, we find that the Hubbard-I approximation with 
practically unrenormalized parameters (t = 0.95 1 and 
U = 8.32 1) fits the QMC spectral functions almost per- 
fectly regarding both the general dispersion and the dis- 
tribution of spectral weight (quadratic deviation per de- 
gree of freedom = 0.05). This is not surprising 
since the Hubbard-I approximation is derived for the 
paramagnetic state thereby neglecting all effects of spin- 
correlations, i.e., the Hubbard-I approximation essen- 
tially describes the interplay between itinerant electrons 
and strong on-site repulsion. This should be a reasonable 
assumption at this high temperature since all relevant 
spin-degrees of freedom are thermally excited. 
Lowering the temperature to T = 1.00 1 and further to 
T = 0.33 1, the data show that the Hubbard-I approxima- 
tion increasingly fails to reproduce the entire spectrum 
and is only able to fit the peaks with maximal spectral 
weight reasonably well. Even then, in order to achieve 
these fits, one already has to renormalize the free param- 
eters strongly. The values we found are t — 1.38 1 and 
U = 5.57 i with = 1.54at T = 1.00 1 and i = 1.40 i and 
U = 3.96 1 with = 0.85 at T = 0.33 <. The peaks that 
are missed by the Hubbard-I approximation are the states 
that form the first ionization/affinity states around mo- 
mentum k = (0, 0)/(7r, tt) on the photoemission/inverse 
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photoemission side and two rather dispersionless bands 
at higher energies of w « ±6.0 The former states were 
previously resolved by Preuss and co-workers All- 
together one can distinguish a total of four 'bands' in 
the single particle spectral density. As will be seen be- 
low, the temperature/doping regime where this 4-band 
structure is seen coincides with the regime where a col- 
lective low-energy mode with momentum (tt, tt) in the 
spin-response exists. Inspite of this, however, we stress 
that the 4-band structure cannot be explained by a back- 
folding of the band structure due to ordering effects since 
the spin-correlation length is < 1.5 lattice spacings at 
temperatures T > 0.33 1. 

For the lowest temperature, T = 0.10 i, the QMC data 
are compared with the results from the AF SDW approx- 
imation. As in the case of the Hubbard-I approximation 
at medium temperatures, the lowermost spectra of fig- 
ure |l| show that the SDW approximation is only able to 
fit the peaks with large spectral weight. Again one has 
to renormalize the free parameters heavily to values of 
i = 1.34 1 and A = 2.29 < with = 0.83. Moreover, as 
was the case for the Hubbard-I approximation at higher 
temperatures, the SDW approximation neither explains 
the states that form the first ionization/affinity states 
around momentum k = (0, 0)/(7r, tt) on the photoemis- 
sion/inverse photoemission side, nor the two dispersion- 
less bands at higher excitation energies which can be seen 
rather clearly in the spectra. 

All in all, the overall distribution of spectral weight is 
roughly reproduced by the Hubbard-I and SDW approx- 
imations as long as one forgets about the 4-band struc- 
ture. In fact it is well known that the integrated pho- 
toemission or inverse photoemission weight (that means 
the electron momentum distribution) at each k-point is 
reproduced quite well by the Hubbard-I approximation 
and the related 2-pole approximation . 
As already mentioned, the emerging of the 4-band struc- 
ture in the photoemission somewhere in between T — 
1.00 i and T = 0.33 i is closely related to a change in 
the spin-response: to illustrate this we consider figure |^, 
which shows the spin-correlation function, XszO^,^) (left 
column), and the charge-correlation function, XccO^,^) 
(right column), for different temperatures. Whereas the 
spin-response is entirely incoherent at T = 1.00 with 
decreasing temperature it can be fitted increasingly well 
by the spin-wave dispersion 



E^^{k) = 2jy 1 - -(cos(fc:,) + cos{ky)y. (19) 

This result is known from previous calculations [Q,^, 
which demonstrated that the two-particle correlation 
functions like the spin-response can be described within 
the SDW approximation even for large values of the in- 
teraction U. The energy scale J directly manifests itself 
in the spin-response since the spin- wave dispersion takes 
the value of E^^^tt.O) = 2 J at momentum k = (tt, 0). 



The fit parameters are J = 0.33 i with x^ = 0-01 at 
T = 0.33t and J = 0.49t with x^ = 0.11 at T = O.lOi. 
The latter is already quite close to the strong coupling 
estimate J = 4<^/J7 = 0.5i. Furthermore, the fig- 
ure shows that with decreasing temperature the spin- 
response concentrates its weight more and more at the 
AF momentum Q — (tt, tt) (as it is the case in AF spin- 
wave theory) and at a characteristic energy uj* . The lat- 
ter decreases with decreasing temperature, i.e., the spin- 
response comes closer and closer to the predictions of 
AF spin- wave theory (^^. The spin-correlation length 
Csz{T) can be derived from a real-space fit of the QMC 
equal-imaginary-times spin-correlation function x(r) to 
the form a ■ |r|'' • e"''''/''"*^-^) thereby incorporating the 
effects of the periodic boundary conditions. While this 
is the best one can do on a finite lattice with periodic 
boundary conditions, the fit will only lead to roughly cor- 
rect values due to the relative small system size of 8 x 8. 
The values obtained for the spin-correlation length then 
are = 0.3 - 0.5, (sz = 1-0 - 1.3, (sz = 1.6 - 1.9, 
Csz = 2.1 - 2.8 and Csz > 8 for T = 1.00 <, T = 0.33 i, 
T = 0.25 1, T = 0.20 1 and T = O.lOt, respectively We 
actually believe that the spin-correlation length reaches 
the system size already at a temperature of T w 0.20 1, 
because at this temperature the fit results in values be- 
tween = 2.1 and Csz = 2.8 (with the exponent b set 
to zero or not) but always with error bars of roughly the 
system size. 

The charge-response Xcc(k, w), on the other hand, is 
rather broad in both momentum k and energy uj for all 
temperatures studied. Furthermore, the charge-response 
is gapped for temperatures below T « 1.00 i and, there- 
fore, can certainly not be responsible for any low-energy 
features of the single particle spectrum. 
It is then quite obvious that at roughly the same tem- 
perature where the two narrow dispersive quasiparticle- 
like bands (that cannot be interpreted within the frame- 
work of the Hubbard-I or SDW approximations) appear 
in the single particle spectrum, the spin-response devel- 
ops a sharp collective low-energy mode. We conclude 
that the underlying mechanism behind the occurrence 
of the 4-band structure consists in dynamical magnetic 
correlation effects, which are beyond the scope of the 
Hubbard-I and SDW approximations. 
In the following, we want to explore the single particle 
spectrum by means of our diagnostic operators. The first 
of them is the shadow operator Ci^a- of equation ( |l5|) which 
will be used to transfer spectral weight from the inner 
parts of the Brillouin zone to the outer ones in case of nor- 
mal photoemission w < /i. As already discussed above, 
this also improves the resolution of the ME method in 
this region, since its resolution strongly depends on the 
spectral weight at a certain position. Nevertheless the 
spectrum of the shadow operator has to exhibit exactly 
the same peak positions as the normal photoemission 
spectrum. 
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Figure ^ shows the angle-resolved spectral function of 
the shadow operator for moderate and low temperatures. 
As expected, the shadow operator has its main spectral 
weight near k = (tt, tt) on the photoemission side and 
near k — (0, 0) on the inverse photoemission side. Fur- 
thermore, it's spectrum supports the existence of a to- 
tal of four bands, because it resolves a group of peaks 
forming dispersionless bands at energies of u; « ±6.0 
a region where the normal photoemission spectrum ex- 
hibits only some weak and smeared-out spectral weight. 
These two dispersionless bands at energies of w « ±6.0 < 
are inconsistent with the dispersions of the Hubbard-I 
and SDW approximations of figure |l|. We will further 
address this topic later in this work. 
Next, we turn in more detail to the temperature depen- 
dence of the photoemission spectrum. Figure ^ shows 
some close-ups of the normal photoemission spectrum 
and of the spectrum of the shadow operator at momen- 
tum k = (tt, tt) and different temperatures. For the nor- 
mal photoemission operator these close-ups show a peak 
at uj fa —1.5 1, which would be consistent with Hubbard-I 
(see Figure ^). In the spectrum of the shadow operator 
this feature is visible as a single resolved peak only at 
the highest temperature, T — 1.00 i whereas for tem- 
peratures down to T = 0.20 i there is only some diffuse 
weight at this position. In the ordinary photoemisson 
spectrum the peak looses spectral weight with decreas- 
ing temperature. It disappears completely at T < 0.20 1 
where the spin-correlation length C,sz {T) reaches the sys- 
tem size (see above). Thus, the temperature T « 0.20 i 
where we lose the Hubbard- I-like peak at w « — 1.5 i and 
k = (tt, tt) coincides quite accurately with the tempera- 
ture where 'effective' long-range order sets in. Moreover, 
we find that as the normal photoemission spectrum loses 
the peak at w ~ —l.5t and k = (7r,7r), the spectrum of 
the shadow operator gains weight at uj « —6.0t. Thus, 
we expect that both features are closely related to the 
temperature development of the spin-correlation length 
Csz(T). We note, however, that the crossover in the shape 
of the dispersion from Hubbard- I-like to SDW-like occurs 
in a quite unexpected way: the topmost band at (tt, tt) 
does not deform into the SDW form in any continuous 
way, but simply 'fades away' and eventually vanishes at 
the transition. 

A further surprising result is the following: at T = O.lOt 
neither the ordinary electron operator nor the shadow 
operator pick up the 'AF umklapp band' correspond- 
ing to the narrow dispersive band seen for example at 
UJ « —3t for k ~ (0,0), i.e. there is no corresponding 
band at a; w —3t and k — {tt,7t). Note that in the 
framework if the SDW-approximation the shadow opera- 
tor must reproduce this umklapp-band at (tt, tt) with the 
same weight as the original band at (0, 0) in the ordi- 
nary photoemission spectrum - see the discussion in the 
first section. That this is not the case shows that even 
at this lowest temperature, a simple SDW-like descrip- 



tion of the band structure is invalid, in that the band 
structure cannot be understood by simple backfolding of 
the spectrum obtained without broken symmetry. As we 
will see in the following the AF SDW state provides only 
the 'background' for the dressing of the photoholes with 
AF spin-excitations which dynamically generate a total 
of four bands. 



III. THE 4-BAND STRUCTURE 

We return to the discrepancy between the 2-band dis- 
persions of the Hubbard-I and SDW approximations and 
the 4-band structure actually observed for example in the 
spectrum of the shadow operator of figure ^ 
In order to generate a '4-band structure' out of the two 
bands of the Hubbard-I and SDW approximations we try 
as a phenomenological ansatz to mix the dispersions of 
the Hubbard-I/SDW approximation with two dispersion- 
less bands at energies of E± = ±3.0 1. In other words, for 
both the photoemission and inverse photoemission spec- 
trum we diagonalize an 'effective' 2x2 Hamilton matrix: 



H± = 



E. 



Hub-I/SDW 

V 



V 

E± 



(20) 



and plot in figure ^ the four bands obtained in this way 
on top of the spectral density obtained from QMC at 
T = 0.33 1 and T — 0.10 1. For comparison the figure also 
shows the original (i.e. unhybridized) Hubbard-I bands 
plus the two phenomenological dispersionless bands at 
energies ol u; — ±3.0 1. The figure shows that the overall 
agreement between the QMC peak positions and the four 



bands generated by diagonalizing H± of equation (20) is 



surprisingly good, particularly so in view of the fact that 
for both, Hubbard-I and SDW approximation, only un- 
renormalized parameters were used. In particular, the 
self- consistently determined value for the SDW gap A of 
3.56 1 was used at T = 0.10 1. The only 'external' param- 
eter in this figure is the mixing matrix element V ^ which 
was set to a value of 1.0 1. 

Thus, we find in contrast to previous works jl4j, that 
the introduction of the dispersionless bands reproduces 
the sinlge-particle gap and the width of the quasiparticle 
band correctly without any renormalization of parame- 
ters. Rather, the narrowing of the quasiparticle band and 
the reduction of the Hubbard gap as compared to the un- 
renormalized parameters is brought about by introduc- 
tion of the dispersionless bands. This naturally raises the 
question as to their physical origin. In the present pa- 
per we restrict ourselves to a more phenomenological and 
'numerics based' approach. A complementary and more 
mathematical discussion is given in Ref. ]ll| ], where an 
equation of motion approach similar to Hubbard's origi- 
nal work is pursued. 

We consider the commutator of the creation operator for 
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hole-like particles, h\ ^ = Cj^(l — nj^s), which annihilates 
a particle only on a singly occupied site, with the kinetic 
energy of the Hubbard model and find : 



1 



(21) 



Keeping only the first term in the square brackets of 
equation (|^) reproduces the Hubbard-I approximation 
pO| . The second term in the square brackets describes 
the dressing of the created hole by a spin-excitation and 
is closely related to the spin- 1/2 string operator of equa- 
tion ( p^l ) . The third term describes in an analogous fash- 
ion the coupling of the hole to a density fluctuations, 
whereas the fourth term describes the coupling to the "q- 
excitation jl^. The two latter types of excitation are 
not important for a large positive U near half-filling, 
(n) = 1.0, and will be neglected. Therefore, the oper- 
ator X]jeAr(j) fS*! + ^S^) is in this case the most 
important correction over the Hubbard-I approximation. 
As already stated, it describes a hole dressed by a spin- 
excitation: this operator not only creates a hole on site 
j but dresses this hole with a spin-excitation on a neigh- 
boring site, which is exactly the idea behind the spin-bag 
or spin-polaron |9| pictures known in the literature. 
Splitting this operator into eigenoperators of Hu- 

qt= fet^f + c,,i5r), 



we find [Di^.Hu] = ^-Dij^^ and [C.ia,Hu] = ~^;Ci,j,<T- 
Assuming moreover that the mobility of these composite 
excitations is determined by the 'heavy' spin-excitation, 
it seems quite resonable to assume that these 'particles' 
are the source of the (more or less) dispersionless bands 
at w ±U/2. 

Finally, the commutation relation (^l|) shows that the 
mixing matrix element between the hk.a and the new 
composite particles should be « t. Based on these rough 
considerations we might thus expect that the two string- 1 
'effective particles' (^ ) are excellent candidates for ex- 
plaining the two dispersionless bands at ±3 1 required to 
upgrade the Hubbard-I or SDW approximation so as to 
match the QMC data. However, so far the above con- 
siderations are pure speculation and in the following we 
will turn to QMC-results to back up this hypothesis by 
numerical evidence. 

Before doing so, however, we want to illustrate the action 
of the 'string-operator' in two extreme cases: an ideal 
Neel state and a resonating valence-bond (RVB) state, 
i.e. a compact singlet covering of the plane [Q (see figure 
PI). In the Neel state, a hole created initially on site i can 



travel one place to a neighboring site j thereby leaving 
behind a misaligned spin on the original site i. Exactly 
this process is described by the second term, ^S^: it 
creates a hole of opposite spin on a neighboring site j and 
flips the spin on the original site i. Therefore, this process 
corresponds to the creation of a string of length 1. In fact 
one might think about more sophisticated diagnostic op- 
erators incorporating the effects of longer-ranged strings 
p6[ . Indeed, Dagotto and Schrieffer and Eder and 
Ohta ||l^ already measured the angle-resolved spectrum 
of a diagnostic operator containing strings with up to 
three lattice sites range by means of exact diagonaliza- 
tions of the t — J model . As already mentioned in the 
introduction, in the QMC method each observable has to 
be expressed in terms of free single-particle Greens func- 
tions on each time slice by the application of Wicks theo- 
rem . This results already in a quite large expression 
for the spin-1/2 string operator of equation (p7| ) contain- 
ing approximately 80 contributions. The implementation 
of even longer-ranged string operators therefore was not 
possible. 

Returning to the spin-1/2 string operator of equation 
( [TtI ) we note that the first term, ^S^ will always an- 
nihilate a Neel state. This reflects the simple fact that 
spin-rotation symmetry is broken in the Neel state. This 
is not the case, however, in the fully rotationally invariant 
RVB state: again, creating a hole on site i and allowing it 
to hop to site j will produce a spin-excitation. However, 
in the case of the RVB state, it produces the superposi- 
tion of two states: in one case the dotted ellipse stands 
for the Sz = component of the triplet, and this state 
(22) would again be created by the term . There is, 

however, also a second state where the dotted ellipse cor- 
responds to the superposition of a singlet and the 5*2 = 
component of the triplet. This second state then would 
be created by the term c^ ^Sj, whereby the relative sign of 
the two terms in the string- 1 operator makes sure that the 
two configurations are always produced with the proper 
phase. In both extreme cases, Neel state and 'singlet 
soup', the string-1 operator thus creates a hole dressed 
with the proper spin-excitation: this scan be a spin wave 
(i.e. a single inverted spin) in the case of a Neel state, or 
a singlet-triplet excitation in the case of an RVB state. 
As a technical remark we still note that the excessive nu- 
merical effort which would have been necessary to com- 
pute spectra for the C'i^a- and Di^^ (which are products of 
5 Fermion operators) has made it impossible to compute 
the spectra of these operators - instead we have been us- 
ing the (Fourier transform of) the operator ^, defined 



in (b_7|). Concerning the difference between this diagnos- 
tic operator and the operator J2jeN{i)i'^j I'^i + S i^T) 
obtained by commuting /i| ^ with the kinetic energy (see 
the second term on the r.h.s. of equation (21)), we note 
that their Fourier transforms differ only by phase factors 
of the form e-«k(Ri-R3)^ Both operators are indeed iden- 
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tical at momentum k = (0, 0) and differ at momentum 
k — (tt, tt) only by a factor of —1 since Ri and Rj are 
next-nearest neiglibor lattice sites. At all other momenta 
both operators have exactly the same peaks but will dif- 
fer somewhat in their spectral weights. 
After these remarks we discuss the angle-resolved spec- 
tral function of the spin- 1/2 string operator 6^. ^ shown 
in figure ^ The spectrum of the spin-1/2 string operator 
indeed picks up those band portions which, according 
to the rough hybridization scenario in figure should 
have strong 'flat-band character', i.e. the part of the 
narrow low energy band between (0,0) and (7r/2,7r/2) 
and between (0, 0) and (tt, 0) at energies w « —3t. We 
note that these are precisely the positions where the two 
quasiparticle-like dispersive narrow bands occurred in the 
normal photoemission spectrum at T = 0.33 1. In addi- 
tion, the band portion at w « —6t for momentum (tt, tt) 
is also enhanced in the string- 1 spectrum. 
This, however, still leaves an important part of the band 
structure unexplained. Namely, the 'AF umklapp band' 
of the narrow quasiparticle band dispersing upward be- 
tween (0,0) and (7r/2,7r/2) at w « —3t still is not seen 
in any of the spectra, not even at the lowest temperature 
studied. On the other hand, in a state with true an- 
tiferromagnetically broken symmetry we know that this 
mirror image must exist due to the backfolding of the 
Brillouin zone. To finally resolve this part, we now intro- 
duce the last diagnostic operator of this work, which we 
call the spin-3/2 string operator: 

^ {2c^^^S]-c,^^S-). (23) 

jeN(t) 

This describes again a composite object of a hole and 
a spin-excitation, but this time the two constituents are 
coupled to the total spin of 3/2. We stress that this oper- 
ator will detect states which can never be seen in an ac- 
tual angle-resolved photoemission spectrum on a singlet 
ground state, because this is forbidden by the angular- 
momentum selection rule. The angle-resolved spectral 
function A(\i,Lu) of the spin-3/2 string operator is plot- 
ted in figure ^ again for T = 0.33 1 (top) and T = 0.10 1 
(bottom). It is then immediately obvious that it is this 
operator which resolves the 'missing piece' of the AF dis- 
persion, i.e. the 'AF mirror image' of the narrow quasi- 
particle band. It should also be noted that the spectrum 
of 6j I is remarkably independent of temperature, i.e. the 
states belonging to this 'spin-3/2 band' persist irrespec- 
tively of whether there is long-range order or not. 
Combining the information obtained so far, suggests the 
following scenario for the crossover between the para- 
magnetic band structure at high temperature and the 
AF band structure at low temperature: in the paramag- 
netic state at high temperatures (such as T = 0.33 1), the 
spin is a good quantum number and the spin-3/2 'band' 
does exist but cannot mix with any spin-1/2 band due 



to spin-conservation. The band of spin-3/2 quasiparticles 
thus plays no role whatsoever in the actual photoemission 
spectrum, which is presumably the reason why the outer 
part of the spectrum is so remarkably 'invisible' in actual 
ARPES experiments, leading to the idea of a 'remnant 
Fermi surface' in the insulator [|l9|. Reaching this state 
by photoemission would only be possible if the photohole 
is created in a thermally admixed state of at least spin 
1. In an infinite system this situation changes discon- 
tinuously at the transition to the true broken symmetry 
state: there the total spin ceases to be a good quantum 
number, and the spin-1/2 band in the interior of the AF 
zone and the spin-3/2 band in the exterior now suddenly 
can mix with each other, thus leading to the familiar 
SDW dispersion. We note that another way to generate 
a coupling between these two bands would be application 
of a magnetic field - this also would break spin-rotation 
invariance and hence enable the hybridization of a spin- 
1/2 and a spin-3/2 band. Based on our results we thus 
believe that a magnetic field could enhance the spectral 
weight of the 'shadow part' of the band structure as seen 
in ARPES. 



IV. DOPING THE HUBBARD MODEL AT HIGH 
TEMPERATURES - RIGID BANDS 

Summarizing our results so far, we may say that the 
Hubbard-I approximation, slightly improved by the intro- 
duction of new quasiparticles corresponding to dressed 
holes provides a very good description of the spectral 
function for the case of half- filling, (n) = 1. In this sec- 
tion we want to proceed to the doped case (n) < 1, which 
is of prime interest for cuprate superconductors. Here, an 
essential drawback of the QMC procedure is that reliable 
QMC simulations for lower temperatures are much more 
difficult or even impossible, since the absence of particle- 
hole symmetry away from half-filling introduces the no- 
torious minus-sign problem into the algorithm. Truely 
low temperatures like T — 0.10 i, which in principle cor- 
respond to the physical temperature range, are therefore 
out of reach. On the other hand in the study of the half- 
filled case we have seen that a major change takes place 
as the spin correlation length reaches the system size, 
whence 'effective long range order' sets in. In the doped 
case the spin correlation length is expected to be short at 
any temperature, whence we may expect that the change 
of A(k, uj) from high to small temperature is more smooth 
than at half-filling. In that sense, even A(k, oj) data for 
the relatively high temperature T — 0.33 i are interest- 
ing to study. Moreover, we can at least try to elucidate 
trends with decreasing temperature and thus construct a 
reasonably plausible scenario. 

At half-filling, we have seen that the 'approximation of 
choice' for the paramagnetic case was the Hubbard-I ap- 
proximation. This naturally poses the question as to how 
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relevant the half-filled case is for the description of the 
doped case, i.e. how much of the Hubbard-I physics re- 
mains valid for finite doping. At half-filling the two 'effec- 
tive particles' dj ^ and cj^., form the two separate Hub- 
bard bands. The effect of doping would now consist in 
the chemical potential cutting progressively into the top 
of the lower Hubbard band, in much the same fashion 
as in a doped band insulator. On the other hand, for 
finite U/t the spectral weight along this band deviates 
from the free-particle value of 1 per momentum and spin 
so that the Fermi surface volume (obtained from the re- 
quirement that the integrated spectral weight up to the 
Fermi energy be equal to the total number of electrons) 
is not in any 'simple' relationship to the number of elec- 
trons - the Luttinger theorem must be violated. This is 
the major reason why the Hubbard-I approximation has 
been dismissed by many authors as being unphysical. 
We now wish to address the question as to what really 
happens if a paramagnetic (i.e. not magnetically or- 
dered) insulator is doped away from half-filling, by QMC 
simulation. We therefore choose T = 0.33t and U — 8.0 1. 
Figure ^ then shows the development of A{k,u!) with 
doping. It is quite obvious from this Figure that initially 
the 2 bands seen at half-filling in the photoemission spec- 
trum (i.e. uj < 0) persist with an essentially unchanged 
dispersion. The chemical potential gradually cuts deeper 
and deeper into the topmost band, forming a hole-like 
Fermi surface centered on (7r,7r), the top of the lower 
Hubbard band. The only deviation from a rather simple 
rigid-band behavior is an additional transfer of spectral 
weight: the part of the topmost band near (tt, tt) gains 
in spectral weight, whereas the band with higher bind- 
ing energy looses weight. In addition, there is a transfer 
of weight from the upper Hubbard band to the inverse 
photoemission part below the Hubbard gap. This effect 
is actually quite well understood |20|. The band struc- 
ture above the Hubbard gap becomes more diffuse upon 
hole doping in that the rather clear two-band structure 
visible near {tt, tt) at half-filling rapidly gives way to one 
broad 'hump' of weight. Apart from the spectral weight 
transfer, however, the band structure on the photoemis- 
sion side is almost unaffected by the hole doping - the 
dispersion of the quasiparticle band becomes somewhat 
wider but does not change appreciably. In that sense we 
see at least qualitatively the behavior predicted by the 
Hubbard I approximation. 

Next, we focus on the Fermi surface volume. Some care 
is necessary here: first, we cannot actually be sure that 
at the high temperature we are using there is still a well- 
defined Fermi surface. Second, the criterion we will be 
using is the crossing of the quasiparticle band through 
the chemical potential. It has to be kept in mind that 
this may be quite misleading, because band portions with 
tiny spectral weight are ignored in this approach (see for 
example Ref. |^l|] for a discussion). When thinking of 
a Fermi surface as the constant energy contour of the 



chemical potential, we have to keep in mind that por- 
tions with low spectral weight may be overlooked. On 
the other hand the fact that a peak with appreciable 
weight crosses from photoemission to inverse photoemis- 
sion at a certain momentum is independent of whether we 
call this a 'Fermi surface' in the usual sense, and should 
be reproduced by any theory which claims to describe 
the system. It therefore has to be kept in mind that in 
the following we are basically studying a 'spectral weight 
Fermi surface', i.e. the locus in k space where an appar- 
ent quasiparticle band with high spectral weight crosses 
the chemical potential. With these caveats in mind. 
Figures |l^ and |ll| show the low-energy peak structure of 
A{k,u!) for all allowed momenta of the 8x8 cluster in 
the irreducible wedge of the Brillouin zone, and for differ- 
ent hole concentrations. In all of these spectra there is a 
pronounced peak, whose position shows a smooth disper- 
sion with momentum. Around (tt, tt) the peak is above 
fi, whereas in the center of the Brillouin zone it is below. 
The locus in fc-space where the peak crosses /x forms a 
closed curve around (tt, tt) and it is obvious from the Fig- 
ure that the 'hole pocket' around (tt, tt) increases very 
rapidly with S. To estimate the Fermi surface volume Vp 
we assign a weight of 1 to momenta k where the peak 
is below fi, 0.5 if the peak is right at /i and if the peak 
is above fi. Our assignments of these weights are given in 
Figures ^ and ^ The fractional Fermi surface volume 
then is Vp — ;^ X^k '"'k, where TV = 64 is the number of 
momenta in the 8x8 cluster. Of course, the assignment of 
the involves a certain degree of arbitrariness. It can 
be seen from Figures ^ and |l^, however, that our Wk 
would in any way tend to underestimate the Fermi sur- 
face volume, so that the obtained Vp data points rather 
have the character of a lower bound to the true Vp ■ Even 
if we take into account some small variations of Vp due to 
different assignments of the weight factors, however, the 
resulting Vp versus S curve never can be made consistent 
with the Luttinger volume, see Figure ^ The deviation 
from the Luttinger volume is quite pronounced at low 
doping. Vp approaches the Luttinger volume for dopings 
~ 20%, but due to our somewhat crude way of determin- 
ing Vp we cannot really decide when precisely the Lut- 
tinger theorem is obeyed. The Hubbard I approximation 
approaches the Luttinger volume for hole concentrations 
of w 50%, i.e. the steepness of the drop of Vp is not 
reproduced quantitatively. The latter is somewhat im- 
proved in the so-called 2-pole approximation jT^,^ . For 
example the Fermi surface given by Beenen and Edwards 
p2t for (n) — 0.94 obviously is very consistent with the 
spectrum in Figure pT| f or (n) — 0.95. 
We return to Figure M and discuss the entire width of 
the spectra, in particular the question of the fate of the 
4-band structure in the doped system. For (n) = 0.95 
the different features that are seen at (n) — 1.0 are 
still rather clearly visible, but for (n) = 0.86 the low 
energy quasiparticle band at A; = (0, 0) starts to dis- 
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appear, and at (n) = 0.80 the dominant 'band' in the 
spectrum between lo = —At and oj = 2t can be fitted 
by a sligtly renormaUzed free-electron band. As we have 
seen above, the Luttinger theorem also is valid in this 
case. This suggests to classify the doping as 'underdopcd' 
for < (n) < 0.85, where the Luttinger theorem is in- 
valid and the 4-band structure known from half-filling 
persists, and 'overdoped' where the Luttinger theorem is 
valid and a renormalized free-electron band can be seen 
in the spectral function. Following the convention for 
cuprate superconductors, we call the doping where the 
crossover between the two regimes occurs the 'optimal' 
doping. 

Next, the four plots of figure |l^ show the spectra at se- 
lected fc-points. The system size is only 6 x 6 in this case 
because this allows for smaller error bars. Closer inspec- 
tion, especially, of the peaks at momentum k = (0, 0) 
on the photoemission side (plot (a)) and at momentum 
k ~ (7r,7r) on the inverse photoemission side (plot (d)) 
confirms, that with increasing hole concentration we are 
losing parts of the 4-band structure seen at half-filling. 
To check the physics of the band structur in more de- 
tail, we again employ our diagnostic operators. Figure 
shows the angle- resolved spectral functions A(k, uj) 
of the spin-1/2 and spin-3/2 string operators, h^-^ and 

&j |. As was the case at half-filling, the spectrum of the 
spin-1/2 string operator highlights exactly those peaks 
that we associate with the dipsersionless 'dressed hole' 
bands in Figure ||. The spectrum of the spin-3/2 string 
operator on the other hand has its peaks with maximal 
spectral weight around momentum k = (tt, tt), indicating 
that also in the doped case there is an 'antiferromagnetic 
mirror image' of the quasiparticle band (which, however, 
consists of spin-3/2 states). Again, coupling of photo- 
holes to thermally excited spin excitations may make 
these states visible in ARPES spectra, thus explaining 
the 'shadow bands' seen in photoemission experiments 
by Aebi et al. [|2^ . Similarly as for half-filling one might 
speculate that a magnetic field, which would break spin 
symmetry and thus allow for a coupling of 'bands' with 
different total spin, would enhance the spectral weight of 
these shadow bands. 

All in all we have ssen that the 'band structure' (4-band 
structure, dispersion of regions of large spectral weight, 
'character' of the bands as measured by the diagnostic 
operators) stays pretty much unchanged as long as we 
are in the underdoped regime. At half-filling the 4-band 
structure is closely related to the sharp low-energy mode 
in the dynamical spin correlation function, which natu- 
rally suggests to study the spin reponse also as a function 
of doping. Figure |l^ shows the spin-correlation func- 
tion, Xs2(k, w) (left column), and the charge-correlation 
function, Xcc(k, w) (right column), for T = 0.33 1 and 
densities (n) = 0.95 (underdoped), {n) = 0.90 (nearly 
optimally doped) and {n) = 0.80 (overdoped). The 



spin-response is sharply confined in both momentum 
k = (tt, tt) and energy lo = uj* only in the underdoped re- 
gion i.e. the regime where we also observe the features as- 
sociated with spin excitations in the single-particle spec- 
tra. As was the case at half-filling for temperatures be- 
low T K, 0.33 1, the spin-response can be fitted by the 
AF spin-wave dispersion (|l^) in the underdoped regime. 
On the other hand, as soon as the system enters the 
overdoped regime the spin-response is no longer sharply 
peaked at momentum k = (tt, tt) and energy ld ^ lo*: it 
broadens in momentum and spreads in energy by an or- 
der of magnitude with the scale changing from J = At'^/U 
to Ekin 8.0 i accompanied by a similar change in the 
bandwidth of the single particle excitations. This re- 
sult is already well known from previous QMC calcula- 
tions ]2^ ] and consistent with similar behaviour in the 
t-J model The charge-response, Xcc(k, w), is always 
broad in both momentum k and energy ui for all densities 
studied. It merely decreases its width from « 12.0 i at 
(n) = 0.95 to « 8.0 1 at (n) = 0.80. 

Although the minus-sign problem of the QMC algorithm 
prevents reliable simulations of large systems at low tem- 
peratures in the doped regime, we nevertheless studied 
the temperature evolution of the angle-resolved spectral 
function A(k, lo) at density (n) — 0.93. This was possible 
due to the relative small system size of 6 x 6, which al- 
leviates the minus-sign problem as compared to 8 x 8 at 
the T — 0.25 1. Figure |l6| shows the results from this 
analysis: the uppermost plot (a) compares the angle- 
resolved spectral functions A(k,uj) at density (n) = 0.93 
for T = 0.50 T = 0.33 i and T = 0.25 i. We stress 
that the simulation at T = 0.25 1 suffers from minus-sign 
problems with a drastically reduced resolution. In the 
center plot (b), the quasiparticle peak weights around 
momentum k = (tt, tt) of the (n) = 0.93 simulation are 
compared with the quasiparticle peak weight at momen- 
tum k = (tt, tt) of a half-filled, (n) = 1.0, simulation for 
different temperatures. At half-filling, the Hubbard-I-like 
quasiparticle peak k — (tt, tt) and lo « —1.5t decreases in 
spectral weight with decreasing T and disappears as the 
spin-correlation length (which increases with decreasing 
temperature) reaches the lattice size (at T « 0.20 1). In 
the underdoped case for density (n) = 0.93 the weights of 
the corresponding peaks around momentum k = (tt, tt) 
located also decrease with decreasing T. Closer inspec- 
tion of this peak (see the inset of the center plot (b)) re- 
veals, that this peak even raises slightly in binding energy 
with decreasing temperature, very similar to the peak in 
the half-filled case. In a real photoemission experiment 
this peak would have dropped below the typical resolu- 
tion of roughly 10% spectral weight ||2^ at a temperature 
of T « 0.25 1. The spin-correlation length (again derived 
by a fit of the equal-imaginary-times spin-correlation 
function to a form a ■ |r|'' • e"''"'/''"'^^^) also shows simi- 
lar behavior in the underdoped and half-filled cases: the 
values for the spin-correlation length are (sz = 0.5 — 0.8, 



10 



(sz = 0.8—1.0 and (sz = 1.0 — 1.3 in the underdoped case 
and Csz = 0.6 - 0.9, (sz = 1-0 - 1.3 and = 1-6 - 1.9 
at half-filling for T = 0.50 1, T = 0.33 1 and T = 0.25 i, 
respectively. The spin-susceptibility (shown in plot (c) 
of figure also behaves very similar, but with changed 
magnitudes in the underdoped and in the half-filled cases. 
These data suggest a similar temperature evolution of the 
band-structure of the Hubbard model in the underdoped 
and half-filled cases driven by the temperature depen- 
dent spin-correlation length (^sziT). Especially, we ex- 
pect the Hubbard-l-like quasiparticle peaks at energies 
of w w 1.0 1 around momentum k = (tt, tt) to vanish with 
decreasing T in the underdoped case as the peaks at en- 
ergies of w w — 1.5t around momentum k = (tt, tt) do in 
the case of half-filling. 

The latter observation suggests a profound change of the 
Fermi surface with temperature: as seen above, it is pre- 
cisely the Hubbard-Tlike band near (tt, tt) which crosses 
the chemical potential and thus forms the Fermi surface 
in the doped case. It is then quite clear that the 'disap- 
pearance' of this band with decreasing temperature must 
affect the Fermi surface in some dramatic way. Studies 
at zero temperature are possible only by means of ex- 
act diagonalization. Analysis of the single particle spec- 
trum shows the same 'rigid-band' behaviour as at high 
temperatures [ p9| and analysis of the momentum distri- 
bution rt(fc) suggest that the doped holes accumu- 
late at the surface of the magnetic zone (i.e. the line 
(f ' f ) ~^ '^)) rather than around (tt, tt). 

V. SUMMARY 

In the present work we have systemactically studied 
the temperature- and doping-dependent dynamics of the 
two-dimensional Hubbard model by finite-temperature 
QMC simulations. Comparing the QMC single parti- 
cle spectral function, the dynamical spin response and 
the spectral functions of suitably chosen diagnostic op- 
erators, different physical regimes could be identified. In 
simplest terms there are two quantities, which basically 
determine the single particle spectrum: the hole concen- 
tration and the the spin-response function, whereby there 
is a certain relationship between the two. 
At half-filling and high temeratures (T > t), the com- 
bined photoemission and and inverse photoemission spec- 
trum A(k, w) displays two dispersive features, the up- 
per and lower Hubbard band, roughly separated by U 
(= 8i, in our work). At these very high temperatures 
the system is in a spin disordered state. We have demon- 
strated here that the well-known Hubbard-I approxima- 
tion gives an excellent description of the single parti- 
cle spectrum in this state, reproducing quantitatively 
both the single-particle dispersion and the distribution 
of spectral weight. This is by no means trivial, since 
the Hubbard-I approximation is dynamically equivalent 



to a simplified effective Hamiltonian, which just contains 
hole-like (hl^) and double occupancy- like particles in a 
simple bi-quadratic form. 

At lower temperatures (T < 0.33 1), the Hubbard-I ap- 
proximation needs to be improved; this is to be expected, 
because it neglects all effects of spin-correlations. In fact, 
the temperature where deviations from Hubbard-I be- 
come strong, coincides fairly well with the transition from 
a spin response Xsz(k, w) which is diffuse both in mo- 
mentum and energy (with a spread of order i), to a more 
'spin- wave-like' response. In this regime Xsz (k,i^) dis- 
plays the characteristic energy scale J = At'^/U, with its 
spectral weight being concentrated at the AF wave vec- 
tor Q = (tt, tt). It should be noted that this 'spin- wave- 
like' regime develops despite the fact that at T = 0.33 1 
the spin-correlations length Csz{T) is still short-ranged 
(< 2 lattice spacings). Only at the lower temperature 
T = 0.10 1 Neel order spreads over the entire QMC block, 
creating an effective (finite-size) Neel state. 
It is well established by previous, in particular also QMC 
work, that in this temperature-regime [T < 0.33 i) new 
spectral features appear. They have often been inter- 
preted as four 'bands', two 'coherent' bands forming the 
topmost valence and the lowest conduction band in the 
insulator plus two 'incoherent' bands, i.e. the remaining 
upper and lower Hubbard band features (see for exam- 
ple |p2[). Our present work not only definitively iden- 
tifies these four bands but also clarifies their physical 
origin and their connection to the spin-excitations. In 
simplest terms the emerging spin waves at lower temper- 
atures provide the excitations that can 'dress' the Hub- 
bard quasiparticles, whence new bands corresponding to 
dressed holes/double occupancies appear in the single 
par ticle spectrum ^(k, w). It has been shown in Ref. 
||ll| that the 4-band structure which appears in A{\i,uj) 
at lower temperatures can be explained in this way, and 
our present numerical check by directly calculating the 
spectra of 'dressed electrons' supports this interpretation. 
This physical picture at half-filling can be extended into 
the underdoped regime. This is most obvious in the single 
particle spectral function, which stays almost unchanged 
in the doped case (i.e. the 4-band structure and the 'char- 
acter' of the bands as measured by the diagnostic oper- 
ators). The main change in fact consists in the chemical 
potential cutting gradually into the (top of the) lower 
Hubbard band, precisely as predicted by the Hubbard-I 
approximation. Contrary to widespread belief the 'Fermi 
surface', if determined by the Fermi surface crossings of 
the dominant band through the chemical potential, does 
not satisfy the Luttinger theorem. Rather, for small hole 
concentrations the Fermi surface volume is considerable 
larger than that for a slightly less than half-filled free- 
electron band. Very similar conclusions have in fact been 
reached by a calculation of the electron momentum dis- 
tribution in the 2D t-J model by Puttika et al. |^l|. Their 
calculation actually was a high temperature series expan- 
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sion plus a Pade extrapolation to lower temperature, and 
it is encouraging that this method gives similar results 
as our QMC results which are performed at relatively 
high temperatures. In its range of applicability, i.e. in 
the absence of strong magnetic correlations and close to 
half-filling, the Hubbard-I approximation thus works re- 
markably well, both at half-filling and in the doped case. 
We stress that this has profound implications for the the- 
oretical treatment of the model: perturbation expansions 
in U or partial and self-consistent resummations thereof, 
may not be expected to give any meaningful results in 
this strong-coupling/low doping regime. 
An interesting question is the possibility to verify our re- 
sults experimentally. As already mentioned above, a scan 
of the temperature development of A{k,uj) shows that 
the part of the quasiparticle band near (tt, tt) (where the 
Fermi surface is located) is loosing weight with decreasing 
temperature. In fact in ARPES experiments on under- 
doped cuprate superconductors the 'hole-pocket' around 
(tt, tt) seen in our simulations (and the expansion of Put- 
tika et al.) is not observed, but rather a small 'Fermi arc' 
near (7r/2,7r/2), terminated by the 'pseudo gap' around 
(tt, 0) Although our simulations do not allow to make 
statements about the truely low temperatures in the ex- 
periments, we believe that this suggests a strong temper- 
ature dependence of the single particle spectrum, with 
the temperature scale being set by the exchange con- 
stant J (which controls the degree of spin disorder) . The 
latter is rather large in copper oxides, so that the tem- 
perature regime studied in our simulations probably can- 
not be accessed experimentally in these materials. We 
note, however, that an ARPES study for the ID mate- 
rial A^ao.96V205 which has a smaller exchange constant, 
has indeed provided evidence for a strong T-dependence 
of A{k^Lo) Clearly, it would be interesting to study 
the Fermi surface evolution in a 2D material with lower 
exchange constant. 

As was the case at half-filling, the dynamical spin- 
response plays an important role: throughout the 
Hubbard-I phase at low doping, the spin response shows 
the sharp low-energy mode at (tTjTt). The simultaneous 
disappearance of the 4-band structure in A{k, oj) and the 
low energy spin response with scale J in the overdoped 
regime then show again the close relationship between the 
two. The dressing of holes by spin excitations apparently 
remains the most important correction over Hubbard-I. 
In the overdoped regime the spin response is spread out 
over an energy range of « St and thus becomes more 
similar to the charge response. The single-particle spec- 
tral function is most consistent with a slightly renormal- 
ized free electron dispersion, and the Luttinger theorem 
appears to be satisfied even at the relatively high tem- 
perature T = 0.33t. This is quite consistent with earlier 
results on the t-J model, which show that for hole concen- 
trations > 25 % the spin and charge response can be ap- 
proximated well by the self-convolution of the single par- 



ticle spectral function [g3|] . This is essentially what is to 
be expected for a system of weakly interacting Fermions, 
so that we conclude that in the overdoped regime we 
enter a new phase which most probably extends to the 
low-concentration limit where the Nagaoka T-matrix ap- 
proximation becomes exact. Finally we note that exact 
diagonalization studies at finite temperatures |Q also 
show some evidence for a 'crossover' between different 
physical regimes at a hole concentration of 15%. 
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FIG. 1. Angle-resolved spectral function ^(k, w) of a 8x8 Hubbard lattice at (n) = 1.0 and U = 8.0 1 for various temperatures. 

The solid linos in the left column represent the (ronormalizcd) Hubbard-I and SDW dispersions. In the right column the spectral 
functions from QMC (solid lines) are compared to the Hubbard-I and SDW- approximations with a Lorentzian lineshape (dashed 
lines). 



FIG. 2. Dynamical spin-correlation function, Xsz(k,a') (left column), and charge-correlation functions Xcc(k, w) (right col- 
umn), of a 8 X 8 Hubbard lattice at (n) = 1.0 and U = 8.0t for T = l.OOt, T = 0.33t and T = O.Wt. The solid lines in the 
lower left column give a spin- wave fit. 



FIG. 3. Angle-resolved spectral function yl(k,a') of the 
shadow-operator „. with 8x8 lattice size at (n) = 1.0 and 
U = 8.0t for s T ='o.33^ and T = O.Wt. 



FIG. 4. Single-particle spectral function of the normal pho- 
toemission (top) and the shadow-operator (bottom) at mo- 
mentum k = (tTjTt) with 8x8 lattice size at (n) = 1.0 
and U = 8.0 1 for temperatures in between T = 1.00 f and 
T = 0.10t 



FIG. 5. Angle-resolved spectral function A(k, ui) of the 8x8 
Hubbard lattice with U = 8.0 1 for T = 0.33 i (top) and 
T = 0.10 1 (bottom) compared to analytic dispersions: mixing 
the Hubbard-I (top, dashed) or SDW (bottom, dashed) bands 
with two dispersionless bands at w = ±3 1 (dashed) results in 
a 4-band structure (solid) in good agreement with the QMC 
peaks. 



FIG. 6. Motion of an aded hole in an Neel state (left) and 
an RVB state (right) produces spin defects of different kind. 
These are described by the spin-1/2 string operator. 



FIG. 7. Angle-resolved spectral function A(k,tj) of the 
spin-1/2 string operator fe^ | in the Hubbard model with 
U = 8.0t for T = 0.33t (top) and T = O.lOt (bottom). The 
lattice size was reduced to 6 x 6 in this case due to the larger 
error bars of the string operator as compared to the normal 
photoemission spectrum. 



FIG. 8. Angle-resolved spectral function A{k.,ui) of the 
spin-3/2 string operator 6^ | in the Hubbard model with 
U = 8.0t for T = 0.33 < (top) and T = 0.10 1 (bottom). The 
lattice size was reduced to 6 x 6 in this case due to the larger 
error bars of the quartet operator as compared to the normal 
photoemission spectrum. 
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FIG. 9. Overview on the doping dependence of the angle-resolved spectral function A(k,aj) of the 8x8 Hubbard lattice at 
T = 0.33* and f/ = 8.0* for densities {n> = 1.00 (half-filled), (n> = 0.95 (underdoped), (n) = 0.86 (roughly optimal doped) 
and (n) = 0.80 (overdoped). 



FIG. 10. Single particle spectral function for all k-points 
of the 8x8 cluster in the irreducible wedge of the Brillouin 
zone. For each k the weight lUk is given. 



FIG. 11. Same as Figure |lO| for lower electron densities. 

FIG. 12. Fermi surface volume as estimated from the single 
particle spectral function, plotted versus the concentration of 
holes in the half- filled band. The dashed line gives the value 
predicted by the Luttinger theorem, Vf ~ ^■ 



FIG. 13. Angle-resolved spectral function j4(k,aj) of a 6 x 6 Hubbard lattice at T = 0.33 1 and U = 8.0 1 for momenta 
k = (0, 0) (plot (a)) and k = (tt, 0) (plot (b)) and k = (tt, tt) (plot (c) and (d)). The spectra for the various densities (n) = 1.00, 
(n) = 0.96, (n) = 0.90, (n) = 0.83 and (n) = 0.66 are shifted by their individual chemical potentials fj, resulting in overlapping 
peaks in the underdoped regime. The labels at the tj-axis refer to the half-filled simulation (solid line). The system size is 
6x6, because the smaller errors in this case lead to more reliable spectra. 



FIG. 14. Angle-resolved spectral function yl(k,a;) of the spin-1/2 string operator 6^ ^ (top) and the spin-3/2 string operator 
(bottom) in the Hubbard model with (/ = 8.0t for T = 0.33 1 and densities (n) = 0.95 (left) and (n) = 0.90 (right). The 
lattice size was 6 x 6 in this case due to the larger error bars of the string operators as compared to the normal photoemission 
spectrum. 



FIG. 15. Dynamical spin-, Xsz (k,tj) (left column), and charge-correlation functions, Xcc(k,Ci;) (right column), of a 8 x 8 
Hubbard lattice at [/ = 8.0* and T = 0.33* for densities (n) = 0.95, (n) = 0.90 and (n) = 0.80. The solid lines in the upper 
plots of the left column give a spin- wave fit. 



FIG. 16. (a) Angle-resolved spectral function j4(k, uo) of the 
6x6 Hubbard lattice with U = 8.0 * and density (n) = 0.93 
versus T, (b) quasiparticle peak weight around momentum 
k = (tt, tt) for (n) = 1.0 (8 x 8) and (n> = 0.93 (6 x 6) versus 
T and (c) magnetic susceptibility Xsz(k = {n,-K),uj — 0) for 
densities (n) = 1.0 (8 x 8) and (n) — 0.93 (6 x 6) versus 
T. Please note the slightly raising (in binding energy) and 
vanishing (in spectral weight) with decreasing temperature of 
the peak located at k = (tt, tt) magnified in the inset of plot 
(b). The system size is only 6 x 6 in the doped case to avoid 
the serious sign problems of the 8x8 system at T — 0.25*. 
But also the 6x6 data at T = 0.25 * suffer from sign problems 
at some momenta, k — (0,0) and k — (7r/3,0). 
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